permeability to estimate
linear barrier permeabilityThe permeability package uses maximum likelihood
estimation with animal biotelemetry data to statistically estimate the
permeability of a linear barrier of interest.
Package functions take an input set of movement tracks with similar sampling rates and quantify crossings of a given linear barrier (Figure 1).
The central function, fitPermeability, uses maximum
likelihood estimation to determine crossing probabilities and to
generate an overall permeability score for the barrier of interest,
together with confidence intervals. Values close to 0 represent a road
that is impermeable and values close to 1 represent a highly permeable
road to animal movement. Values greater than 1 are also possible, with a
“hyper-permeable” barrier.
Additional functions can be used to a) estimate permeability of the barrier as a function of covariate(s), b) compare models, and c) summarize and visualize model results.
This vignette demonstrates use of package functions with both simulated and “real” animal movement data.
To load the package:
Given a set of consecutive animal locations (\(Z\)), “permeability” (\(\kappa\)) is modeled as the influence of a linear element (e.g., a road or barrier, \(B\)) on the animal’s redistribution kernel, aka the probability distribution of finding an animal one-time unit (or “step”) after a given observation. A barrier with a permeability of 0 results in a zero probability density on the other side, whereas a semi-permeable barrier results in a proportion of the distribution, scaled relative to the null distribution, occurring on the other side of the barrier. The amount that the redistribution kernel on the other side of the barrier is scaled relative to the null distribution can be captured with the parameter, \(\kappa\), with values ranging from 0 to 1, or even > 1 in the case of a “hyper-permeable” barrier.
The probability at time t that an animal crosses a given barrier is given by:
\[c_t(Z_t, B, \kappa) = 1 - s_0(Z_t, B)^\kappa\]
where \(s_0(\cdot) = 1 - c_0(\cdot)\) is the null probability that the animal “stays” (does not cross barrier \(B\) from location \(Z_t\)). This formulation has the desired properties: when \(\kappa = 0\), the probability of crossing is 0, when \(\kappa = 1\), the probability of crossing is the null probability \(c_0 = 1-s_0\). As \(\kappa\) approaches infinity, the probability of crossing approaches 1.
Thus, given a set of \(k\) null steps \(\vec{S_i}\) where \(i \in \{1,2,...,k\}\) we take each location of interest \(Z_t\), add the steps \(S_i\), and simply count the number of relocations on either side of the barrier \(B\).
Once this null set is obtained, the number of null locations that crossed and didn’t cross the boundary can be counted (\(n_{c,t}\) and \(n_{s,t}\), respectively) and the estimate of the null probability of staying \(s_0\) found as:
\[\widehat{s_{0,t}} = \frac{n_{s,t}}{n_{c,t} + n_{s,t}}\]
The movement and barrier data can thus be reduced to a fundamental set of interacting observations, stored within a data frame that we define as the “crossing table”, where each row represents a single potential crossing event. At a minimum, the table contains a vector of actual crossing events \(c_i = \{0,1\}\), a column of estimated null crossing probabilities \(\widehat{c_{0,i}}\), and a column of estimated null “staying” probabilities \(\widehat{s_{0,i}} = 1- c_{0,1}\). The crossing table can further be annotated with any number of covariates, as detailed below.
After obtaining the crossing table, the log-likelihood function of \(\kappa\) given movement locations \(Z\) and barrier \(B\), is:
\[{\cal L}(\kappa | Z,B) = \prod_{i=1}^n \left(I(c_i) (1 - \widehat{s_{0,i}}^\kappa) + (1 - I(c_i)) \widehat{s_{0,i}}^\kappa \right)\]
where \(I(c_i)\) is the indicator function of having crossed at time i, i.e. the first term represents the rows of the crossing table (movement steps) where the animal actually crossed the barrier and the second term represents the movement steps where the animal did not cross the barrier.
From this, the (log)-likelihood is readily maximized to estimate \(\kappa\) with approximate confidence intervals (via the Hessian of the likelihood at the maximum likelihood estimator).
The permeability parameter \(\kappa\) can also be estimated as an exponential function of covariate(s) \(X\):
\[\kappa = \exp(\beta X)\]
Covariates may represent a feature of either the barrier (e.g., a static variable, such as variable structure with different barrier segments, or a dynamic variable, such as vehicle traffic or kilometer along the highway) or the animal movement (e.g., a static variable, such as animal sex, or a dynamic variable, such as season). Multiple covariates can be included in a multivariate model formula - however, users should be cautious of parsimony, over-fitting, and sample size. When fitting categorical data, users should ensure that there is a large enough sample size of steps proximate to the barrier to estimate \(\kappa\) for different category levels.
The choices made in creating null steps are important, as the distribution of step lengths sampled to create the null set of steps can impact the null probability of crossing the barrier and therefore, the estimation of \(\kappa\).
Null steps should accurately reflect crossing behavior, such that each null step has the “potential to cross”. This may be best accomplished by simply sampling from an empirical distribution of step lengths for steps that cross and from a distribution of turning angles observed for any steps proximate to the barrier.
This approach accounts for crossing-specific differences in step lengths, commonly observed for animals interacting with anthropogenic features, e.g., animals “bolting” across a road with larger step lengths to reduce risk of harm. Notably, for datasets with a low number of crossings (\(< 30\)), there may not be enough steps that cross to create an empirical distribution of step lengths to sample from and a simulated distribution may be required (similar to methods used for integrated Step Selection Analyses).
Use of the permeability package assumes that:
sf package)and
sf LINESTRING
geometry type or a matrix of barrier segment locations) that the
movement data interacts with.We demonstrate usage of the package and its associated functions with both simulated and “real” movement data below.
Using the simulatePermeability function, a hypothetical
movement track (as a correlated random walk, with a length \(N\), turning angles \(\theta\), and step lengths \(L\)) can be simulated with an example
barrier (composed of a set of linear segments), with the barrier having
a pre-defined \(\kappa\).
The user can mainly define:
Shape and scale parameters (step.shape,
step.scale) for simulated track step lengths (Weibull
distribution, see rweibull function)
Rho or concentration parameter (theta.rho) for track
turning angles (wrapped Cauchy distribution, see
rwrappedcauchy function from the circular
package)
Track length (N.steps)
Bounding box x and y limits for the simulation
(xmax, ymax)
Number of segments in the barrier
(n.segments)
Desired permeability of the barrier (kappa)
If the plot.track argument is set to TRUE, a plot of the
simulated track and barrier will be returned, with steps that crossed
the barrier shown in purple and steps that “bounced” off the barrier
shown in orange:
sim_track <- simulatePermeability(xmax = 30, ymax = 30,
N.steps = 500, theta.rho = 0.5,
step.shape = 5, step.scale = 10,
n.segments = 10,
kappa = 0.2,
plot.track = TRUE) The output of this simulation is a permdata list object,
which contains two elements:
track - data frame containing:
Start and end locations of each movement step
(Z.start,Z.end)
Time difference between consecutive locations
(D.time)
Simulated timestamps for each location
(Time)
Movement step as a complex location (Step, see
as.complex())
Step length a.k.a the Euclidean distance between consecutive
locations (L)
Absolute and turning angles for locations (Phi and
Theta)
A unique identifier for each step (Step.ID)
Whether the step was within the maximum step length of the
barrier (In.buffer) - only steps proximate to the barrier
(where null crossings are possible) are needed to estimate \(\kappa\) - for simulated data from
simulatePermeability, all steps are retained (all simulated
steps have In.buffer==TRUE)
Distance to barrier for the start location in each step
(Dist_toBarrier)
Whether the step crossed the barrier or not
(Crossed) - NA values for steps where
In.buffer==FALSE
Simulated unique track identifier (ID)
And:
barrier - data frame containing:
Complex locations (see as.complex()) for the start
and end of each barrier segment (Z1,
Z2)
Identifier for each line (in sf terms, LINESTRING)
in the barrier (line_id) - important for barriers with
multiple lines (e.g., a grid-like or network-like barrier) - will be set
to 1 if there is only 1 unique line
Identifier for each barrier segment
(barrier.id)
“True” permeability \(\kappa\)
values (kappa) for each barrier segment, based on your
input parameters
If a covariate model is being used, optionally covar
values, one for each unique barrier segment (see “Covariate Model”
example below)
str(sim_track)
> List of 2
> $ track :'data.frame': 499 obs. of 13 variables:
> ..$ Z.start : cplx [1:499] 28.8-8.7i 21.5-11.4i 27.5-22.2i ...
> ..$ Z.end : cplx [1:499] 21.5-11.4i 27.5-22.2i 29.8-13.2i ...
> ..$ D.time : num [1:499] 1 1 1 1 1 1 1 1 1 1 ...
> ..$ Time : POSIXct[1:499], format: "2025-02-25 00:00:00" "2025-02-25 01:00:00" ...
> ..$ Step : cplx [1:499] -7.3-2.69i 6.02-10.8i 2.33+8.99i ...
> ..$ L : num [1:499] 7.78 12.37 9.29 8.67 10.01 ...
> ..$ Phi : num [1:499] -2.79 -1.06 1.32 2.67 2.53 ...
> ..$ Theta : num [1:499] NA 1.726 2.38 1.348 -0.139 ...
> ..$ Step.ID : num [1:499] 1 2 3 4 5 6 7 8 9 10 ...
> ..$ In.buffer : logi [1:499] FALSE FALSE FALSE FALSE FALSE FALSE ...
> ..$ Dist_toBarrier: num [1:499] 28.9 21.5 27.5 29.8 22.3 ...
> ..$ Crossed : logi [1:499] FALSE FALSE FALSE FALSE FALSE FALSE ...
> ..$ ID : Factor w/ 1 level "SimID_1": 1 1 1 1 1 1 1 1 1 1 ...
> $ barrier:'data.frame': 10 obs. of 5 variables:
> ..$ Z1 : cplx [1:10] -0.2-30i 0-24i -0.2-18i ...
> ..$ Z2 : cplx [1:10] 0-24i -0.2-18i 0-12i ...
> ..$ line_id : num [1:10] 1 1 1 1 1 1 1 1 1 1
> ..$ kappa : num [1:10] 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2 0.2
> ..$ barrier.id: int [1:10] 1 2 3 4 5 6 7 8 9 10
> - attr(*, "class")= chr "permdata"
> - attr(*, "buffer.size")= num 5.37
> - attr(*, "sample.rate")= num 1Inmportantly, both track and barrier are in
a segmented step format, i.e. each row represents a segment
(i.e., a movement step or a barrier segment), and it is in the
interaction of the movement segments and the barrier segments that the
permeability is analyzed.
head(sim_track$track)
> Z.start Z.end D.time Time
> 1 28.76543- 8.68803i 21.46299-11.37853i 1 2025-02-25 00:00:00
> 2 21.46299-11.37853i 27.48339-22.18237i 1 2025-02-25 01:00:00
> 3 27.48339-22.18237i 29.80944-13.19165i 1 2025-02-25 02:00:00
> 4 29.80944-13.19165i 22.10698- 9.22025i 1 2025-02-25 03:00:00
> 5 22.10698- 9.22025i 13.93393- 3.44549i 1 2025-02-25 04:00:00
> 6 13.93393- 3.44549i 13.14723+ 8.57085i 1 2025-02-25 05:00:00
> Step L Phi Theta Step.ID In.buffer
> 1 -7.302441- 2.690500i 7.782315 -2.788587 NA 1 FALSE
> 2 6.020405-10.803844i 12.368037 -1.062407 1.7261804 2 FALSE
> 3 2.326045+ 8.990724i 9.286744 1.317631 2.3800378 3 FALSE
> 4 -7.702458+ 3.971395i 8.666016 2.665542 1.3479110 4 FALSE
> 5 -8.173053+ 5.774759i 10.007329 2.526477 -0.1390653 5 FALSE
> 6 -0.786696+12.016338i 12.042062 1.636172 -0.8903051 6 FALSE
> Dist_toBarrier Crossed ID
> 1 28.91244 FALSE SimID_1
> 2 21.46246 FALSE SimID_1
> 3 27.49371 FALSE SimID_1
> 4 29.82373 FALSE SimID_1
> 5 22.27161 FALSE SimID_1
> 6 14.18758 FALSE SimID_1head(sim_track$barrier)
> Z1 Z2 line_id kappa barrier.id
> 1 -0.20501621-30i 0.04983172-24i 1 0.2 1
> 2 0.04983172-24i -0.19873887-18i 1 0.2 2
> 3 -0.19873887-18i 0.00952611-12i 1 0.2 3
> 4 0.00952611-12i -0.02178609- 6i 1 0.2 4
> 5 -0.02178609- 6i -0.04682003+ 0i 1 0.2 5
> 6 -0.04682003+ 0i -0.15139704+ 6i 1 0.2 6In addition to a correlated random walk (type = "crw" in
simulatePermeability), users can also specify a biased
correlated random walk (type = "bcrw"), including an
argument for the target.location (as a matrix of X/Y
coordinates) if they want to set the target location themselves
(otherwise the package will simulate a random target location).
sim_track <- simulatePermeability(xmax = 50, ymax = 50,
N.steps = 100, theta.rho = 0.5,
step.shape = 5, step.scale = 10,
n.segments = 10,
kappa = 0.2,
type = "bcrw",
target.location = target,
plot.track = TRUE)The user can also supply a simulated barrier of their own to the
simulatePermeability function.
For example, a “grid-like” barrier could be simulated:
xmax <- 50
ymax <- 50
n.segments <- 10
barrier_sep <- (1/2)*xmax # spacing between linear features
y <- ymax-barrier_sep # allow spacing around barrier lines
barrier1 <- cbind(X = (xmax-barrier_sep),
Y = y * seq(-1,1,length = n.segments+1),
line_id = 1)
barrier2 <- cbind(X = (xmax-barrier_sep*2),
Y = y * seq(-1,1,length = n.segments+1),
line_id = 2)
barrier3 <- cbind(X = (xmax-barrier_sep*3),
Y = y * seq(-1,1,length = n.segments+1),
line_id = 3)plot(barrier1, type= "l", ylim = c(-ymax,ymax), xlim = c(-xmax,xmax))
lines(barrier2)
lines(barrier3)
grid.barrier <- rbind(barrier1, barrier2, barrier3)
head(grid.barrier)
> X Y line_id
> [1,] 25 -25 1
> [2,] 25 -20 1
> [3,] 25 -15 1
> [4,] 25 -10 1
> [5,] 25 -5 1
> [6,] 25 0 1The simulated barrier can be used instead of generating a random
barrier, by providing it in the barrier argument -
n.segments = NULL, as n.segments is now equal
to the number of segments present in the user-inputted barrier
object.
Importantly, the simulatePermeability function assumes
that for a barrier with multiple lines that there is a equal number of
segments in each line.
sim_track <- simulatePermeability(xmax = 50, ymax = 50,
N.steps = 500, theta.rho = 0.5,
step.shape = 5, step.scale = 10,
kappa = 0.2,
barrier = grid.barrier,
plot.track = TRUE)
The crossing table for the simulated data can be constructed using
the buildCrossingTable function, which takes a
permdata object (containing simulated movement
track and barrier objects) and optionally, the
number of null steps to sample around each observed step (default
n.null = 60).
The resulting object, similar to the track object in
permdata, is a “step-wise” dataframe with rows for every
observed step in the simulated data that had null crossings with the
barrier (e.g., had the “potential to cross”) and with new columns for
the:
Number of null steps to cross/not cross the barrier
(null.crossed, null.stayed)
Number of actual steps in the data to cross the barrier
(crossed)
Barrier id for the barrier segment where the majority of possible
steps crossed (barrier.id - corresponds to the same
barrier.id in permdata)
If the actual step crossed, the barrier segment where it crossed
(barrier.crossed) - barrier.crossed = NA if
the actual step did not cross the barrier
str(sim_cT)
> 'data.frame': 46 obs. of 10 variables:
> $ ID : Factor w/ 1 level "SimID_1": 1 1 1 1 1 1 1 1 1 1 ...
> $ Step.ID : num 4 5 10 12 14 16 113 114 115 120 ...
> $ barrier.id : num 19 18 28 29 25 23 25 27 28 29 ...
> $ Z.start : cplx -3+19i -3.1+10.5i -24.5+13.3i ...
> $ Z.end : cplx -3.08+10.47i -4.79+7.51i -14.67+20.87i ...
> $ Time : POSIXct, format: "2025-02-25 03:00:00" "2025-02-25 04:00:00" ...
> $ null.crossed : int 26 20 27 28 19 12 24 24 25 21 ...
> $ null.stayed : int 34 40 33 32 41 48 36 36 35 39 ...
> $ crossed : num 0 0 0 1 0 0 0 1 0 0 ...
> $ barrier.crossed: num NA NA NA 29 NA NA NA 27 NA NA ...To estimate \(\kappa\) using the
compiled crossing table, the fitPermeability function can
be used.
If no equation is supplied, a null model is fit:
null_model <- fitPermeability(data = sim_cT)
null_model
>
> Estimate of kappa (with 95% confidence intervals):
> kappa.hat ci.low ci.high
> 1 0.2446779 0.1017568 0.5883373The print, summary, and
predict functions can be used to examine the output,
including the Log-Likelihood, AIC, and \(\kappa\) estimate with 95% confidence
interval.
summary(null_model)
> $coef
> estimate se z_score p_value
> 1 -1.407813 0.4476397 -3.144968 0.001661048
>
> $logLik
> [1] -14.43865
>
> $AIC
> [1] 30.87731To potentially get more robust estimates (with tighter standard errors), a larger sample size (e.g., 50 tracks) can be used.
The simulatePermeability function has an additional
argument available for n.tracks, with the default being 1.
The function can also be run with parallel = TRUE (using
the future.apply package) to speed up calculations.
Permeability can also be estimated as a function of covariate(s).
Covariates can be specific to the barrier (e.g., one unique covariate
value per barrier segment, e.g. unique to each barrier.id),
to the movement track (e.g., one unique covariate value for each
Step.ID), or the interaction of the barrier and movement
track (e.g., specific to crossing table rows). We demonstrate covariate
values specific to the barrier here and for the movement track and
interaction of the barrier and movement track in the boreal caribou
“real data” examples below.
The simulateTrack function will simulate covariate
values for you (if no barrier object is specified) or, if a
user-provided barrier object is input, will take an
argument for covars as a vector of numeric values
corresponding to the desired covariate values for each barrier segment
(barrier.id).
Additionally, instead of supplying a kappa value, we
supply the coefficients defining the (log-linear) relationship between
\(\kappa\) and our covariate, with
beta storing the desired intercept (\(\beta_0\)) and slope (e.g., \(\beta_1\)) values, respectively.
The plot shows barrier segments with low covariate values in black and barrier segments with high covariate values in green:
sim_track <- simulatePermeability(xmax = 50, ymax = 50,
N.steps = 500, theta.rho = 0.5,
step.shape = 5, step.scale = 10,
n.segments = 15,
beta = c(-2,1),
plot.track = TRUE)
Note that the barrier object within the
permdata output now also has additional elements for the
covariate values and their associated “true” \(\kappa\) values:
head(sim_track$barrier[,c("covar","kappa")])
> covar kappa
> 1 2.500740 0.03303250
> 2 2.645828 0.03383445
> 3 2.761895 0.03449000
> 4 5.480341 0.05406118
> 5 6.310057 0.06200998
> 6 9.729632 0.10914329The crossing table is built the same as for the null model example:
After building, the crossing table needs to be annotated with the covariate values. This allows flexibility for covariates, based on whether the covariates are associated with the barrier itself or some other non-barrier (e.g., movement step) specific covariate.
For our simulated example, we have one covariate value for each
barrier segment (stored in the barrier object in
sim_track), so we can use the barrier.id
column in the crossing table to annotate our covariate values to our
crossing table.
sim_cT <- merge(sim_cT, sim_track$barrier[,c("barrier.id","covar")],
by = "barrier.id", all.x = TRUE)
str(sim_cT)
> 'data.frame': 19 obs. of 11 variables:
> $ barrier.id : num 1 2 2 4 4 6 7 7 8 8 ...
> $ ID : Factor w/ 1 level "SimID_1": 1 1 1 1 1 1 1 1 1 1 ...
> $ Step.ID : num 392 363 367 379 384 386 137 109 108 216 ...
> $ Z.start : cplx -2.1-42i -1.6-38.3i -3.7-49i ...
> $ Z.end : cplx -8.9-34i -13-37.3i -3.8-44.8i ...
> $ Time : POSIXct, format: "2025-03-13 08:00:00" "2025-03-12 03:00:00" ...
> $ null.crossed : int 16 29 5 40 12 37 39 34 20 11 ...
> $ null.stayed : int 44 31 55 20 48 23 21 26 40 49 ...
> $ crossed : num 0 0 0 0 0 0 0 0 0 1 ...
> $ barrier.crossed: num NA NA NA NA NA NA NA NA NA 7 ...
> $ covar : num 2.5 2.65 2.65 5.48 5.48 ...To estimate \(\kappa\) as a function
of covar, we supply our covariate column from our crossing
table to the formula:
The estimated coefficient values can be extracted with the
coef.table element:
covar_model$coef.table
> estimate se z_score p_value
> (Intercept) -1.581248 0.7986390 -1.979928 0.04771156
> covar 1.196793 0.8500841 1.407853 0.15917465Estimates may be improved by simulating multiple tracks.
The predict function can be used to predict \(\kappa\) values for each barrier segment
and unique covar value.
Be sure to specify se.fit=TRUE to get upper and lower
(95%) confidence intervals on \(\kappa\) estimates with the
predict function.
predict_kappa <- predict(covar_model, se.fit = TRUE)
predict_kappa$covar <- covar_model$data$covar
predict_kappa <- predict_kappa |>
dplyr::arrange(kappa.hat) # sort by increasing kappa values
head(predict_kappa)
> kappa.hat ci.low ci.high covar
> 1 0.02550785 0.0004041883 1.6097707 2.500740
> 2 0.02651591 0.0004425351 1.5887852 2.645828
> 3 0.02651591 0.0004425351 1.5887852 2.645828
> 4 0.05654019 0.0025631688 1.2472034 5.480341
> 5 0.05654019 0.0025631688 1.2472034 5.480341
> 6 0.17593241 0.0314121231 0.9853588 9.729632A plot can be made of the predicted \(\kappa\) values vs the covariate values for our simulated tracks, with 95% confidence intervals.
plot(predict_kappa$covar, predict_kappa$kappa.hat, type = "l",
xlab = "Covariate Values", ylab = "Estimated Kappa Values", col ="orange", lwd = 2, ylim = c(0, 2))
lines(predict_kappa$covar, predict_kappa$ci.high, col = "grey")
lines(predict_kappa$covar, predict_kappa$ci.low, col = "grey")Additionally, if the user wishes to supply their own barrier with
covariate values for each segment, this can be done with the
covar argument.
Boreal caribou (Rangifer tarandus) are non-migratory caribou inhabiting boreal forests in Canada, including areas proximate to human development and major road highways in the Northwest Territories (NWT).
In the Northwest Territories, boreal caribou are known to utilize habitat proximate to major highways such as Highway 1 (Figure 3) but may avoid crossing when vehicle traffic is high.
This example will use GPS tracking data from 10 boreal caribou that encounter Highway 1 West, one of the major NWT highways.
Example data, including a dataframe of caribou locations, Highway 1
West sf object, and daily traffic volume data for Highway 1
West can be loaded using the data function.
The permeability package requires movement data to have
columns for:
ID (factor or character
structure)
Time (POSIXct, date and time)
X and Y spatial coordinates (X and Y,
units in meters)
Any covariate(s) specific to movement steps (can be numeric or
categorical - here, traffic_volume and
Season)
str(boreal_caribou)
> 'data.frame': 9522 obs. of 5 variables:
> $ ID : Factor w/ 10 levels "ID_1","ID_3",..: 1 4 4 1 4 1 1 4 4 1 ...
> $ Time : POSIXct, format: "2015-02-20 17:37:00" "2015-02-20 23:24:00" ...
> $ X : num 330632 409718 427748 327542 427499 ...
> $ Y : num 6787977 6729661 6736940 6784926 6738101 ...
> $ Season: chr "winter" "winter" "winter" "winter" ...For this example, we will focus on caribou interactions with the western portion of Highway 1.
Barrier data can either be a spatial LINESTRING object (class
sf) or a matrix of X and Y coordinates. It should have the
same coordinate reference system as the movement data.
str(hwy1_west)
> Classes 'sf' and 'data.frame': 29 obs. of 1 variable:
> $ geometry:sfc_LINESTRING of length 29; first list element: 'XY' num [1:568, 1:2] 472462 472177 472153 472135 472106 ...
> - attr(*, "sf_column")= chr "geometry"
> - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..:
> ..- attr(*, "names")= chr(0)X and Y coordinates can be extracted from the sf object
and converted to a matrix using the st_coordinates
function.
hwy1_west_xy <- st_coordinates(hwy1_west)
head(hwy1_west_xy)
> X Y L1
> [1,] 472461.7 6771187 1
> [2,] 472176.9 6771241 1
> [3,] 472152.6 6771245 1
> [4,] 472134.7 6771248 1
> [5,] 472106.4 6771251 1
> [6,] 472086.5 6771252 1The output includes a third column L1 for the id of the
line. If the barrier object has multiple lines (e.g., a grid or network
format), this column is especially useful for keeping track of distinct
geometric lines. However, note that even features that appear to be one
line are often made of multiple lines:
unique(hwy1_west_xy[,"L1"])
> [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
> [26] 26 27 28 29If barrier data is being provided as a non-sf object
(matrix of barrier segments), permeability functions
require the barrier data to contain X and Y
coordinates and the third column to be named line_id, even
if there is only 1 line in the object.
Alternatively, users can simply provide their barrier data as an
sf LINESTRING object to permeability
functions.
Users can also annotate these coordinates with covariate values
specific to each segment (see simulation examples above). Just keep in
mind that for each line in your barrier coordinates matrix (e.g., each
unique line_id), there are \(n\) rows but \(1
- n\) segments.
Additionally, if your LINESTRING object is continuous in space but
made of multiple lines (e.g., the hwy1_west example here,
which has 29 lines within a continuous feature), you will need to
account for the fact that the end of each line has the same spatial
point as the beginning of the next line (resulting in duplicates, unless
you remove them).
hwy1_west_xy[566:573,]
> X Y line_id
> [1,] 394810.9 6777956 1
> [2,] 394765.2 6777979 1
> [3,] 394762.2 6777981 1
> [4,] 394762.2 6777981 2
> [5,] 394714.5 6778006 2
> [6,] 394714.5 6778006 3
> [7,] 394066.9 6778348 3
> [8,] 394015.8 6778374 3Examples are shown for this with the covariate model examples below.
It’s helpful to check using a basic plot of the barrier coordinates
(as complex locations, see as.complex) that the barrier
spatial object looks as it should (occasionally, imported shapefiles
from ArcGIS can have unordered segments).
Note that if you have multiple lines in your barrier (e.g., multiple
unique line_id values), each line should be plotted
separately to ensure correct ordering of segments.
The permeability package requires that movement data
have a (mostly) regular fix rate. The fix rate should be checked prior
to utilizing permeability functions.
time_steps <- boreal_caribou |>
group_by(ID) |>
arrange(Time) |>
mutate(dT = c(NA, as.numeric(round(difftime(Time[-1], Time[-length(Time)], units = "hours"))))) |>
ungroup() |>
data.frame() |>
select(dT)The example data is primarily composed of fixes every 8 hours.
Some irregularity is tolerable, since steps not equal to the fix rate
(with a tolerance that can be specified with
sample.rate.tolerance in the
prepPermeabilitymfunction - default 1 hour) will simply not
be used in permeability functions - this ensures that the user does not
need to use interpolation (which may result in “fake” regular
locations).
The user should ensure that their data is clean (processed to remove
any tag errors or other erroneous locations) and mostly regular. We
recommend the resample function from the amt
package to aggregate any steps with short fix rates less than the
desired sample.rate.
Note that if there are many irregular timesteps or large
gaps between locations in your data, the user should regularize their
data before using the permeability package functions, as
otherwise there will not be enough movement steps to estimate \(\kappa\).
After the tracking data has been cleaned, the
prepPermeability function is used to further format
movement data for permeability package functions, creating
a permdata object similar to the simulation examples.
Main arguments to set in this function are the
sample.rate (and optionally,
sample.rate.tolerance, which defaults to 1 hour around the
sample rate) and optionally, the buffer around the barrier
(buffer.m). If buffer.m = NULL (the default),
prepPermeability will set the buffer to the maximum
distance from the barrier observed for steps that cross. The user can
also supply their own buffer.m, but should only do so if
this distance justifiably will include all steps with the “potential to
cross” the barrier (all steps not in the buffer will not be used to
build the crossing table).
Setting plot.track=TRUE will also optionally output a
plot of each track, including the barrier in red if it’s proximate. As
example for one individual:
bor_format <- prepPermeability(move.df = subset(boreal_caribou, ID == "ID_3"),
barrier = hwy1_west_xy,
sample.rate = 8,
plot.track = TRUE)Covariates can also be annotated to the permdata object.
Covariates specific to animal movement steps can be provided with the
move.covar.names and move.covar.values
arguments. Barrier segment specific covariates can be specified using
the barrier.covar.values argument, which should contain a
dataframe of covariate values for every barrier segment (e.g., rows
correspond to each barrier segment or each unique
barrier.id, AFTER accounting for unique lines, e.g., any
duplicate segments).
For movement step specific covariates, the
prepPermeability function takes a vector with the name(s)
of the covariate column(s) in the input movement data
(move.covar.names), as well as an argument for how to
annotate steps with covariate data (move.covar.values,
options being to use covariate values from the start
location in the movement step, the end location in the
step, or for numeric covariates, optionally the mean or
difference of the start and end location covariate values -
defaults start).
For our example boreal caribou data, locations have been annotated with 5 seasons (pre-calving Apr, calving May - June, summer Jul - Aug, fall Sep - Nov, and winter Dec - Apr), a movement-step specific covariate.
We also create a barrier-specific covariate, km,
representing kilometer markers for every barrier segment along the
highway, from the complex locations derived from the barrier X/Y
coordinate data and accounting for unique lines.
z <- hwy1_west_xy[,1] + 1i*hwy1_west_xy[,2]
barrier_covar <- cbind(hwy1_west_xy, z) |>
data.frame() |>
group_by(line_id) |>
reframe(km = Mod(diff(z))) |>
select(km) |>
mutate(km = cumsum(km)/1000) |>
data.frame()To speed up run time, this function can also be run in parallel,
using thefuture.apply package and
parallel = TRUE):
bor_format <- prepPermeability(move.df = boreal_caribou,
move.covar.names = c("Season"),
barrier.covar.values = barrier_covar,
barrier = hwy1_west,
sample.rate = 8)The output permdata object is a list with dataframes
track and barrier.
track contains the movement data, with rows for every
step and including columns for:
Start and end locations of the step (Z.start,
Z.end, as complex locations, see
?as.complex())
Time length of the step (in hours, D.time)
Date/time for the first location in the step (Time,
POSIXct format)
Various movement metrics (step length L in meters,
absolute and turning angles in radians, Phi and
Theta)
Whether the step was within the buffer of the barrier
(In.buffer)
Unique identifier for the step in each individual track
(Step.ID)
Whether the step crossed the barrier
(Crossed)
Any step-specific covariate values (Season)
str(bor_format$track)
> 'data.frame': 9512 obs. of 14 variables:
> $ Z.start : cplx 330632+6787977i 327542+6784926i 327254+6784985i ...
> $ Z.end : cplx 327542+6784926i 327254+6784985i 328011+6787781i ...
> $ D.time : num 14 8 8 8 8 8 8 8 8 8 ...
> $ ID : Factor w/ 10 levels "ID_1","ID_11",..: 1 1 1 1 1 1 1 1 1 1 ...
> $ Time : POSIXct, format: "2015-02-20 17:37:00" "2015-02-21 08:03:00" ...
> $ Step : cplx -3090-3051i -288+59i 757+2797i ...
> $ L : num 4342.6 294.1 2897.1 34.1 39.3 ...
> $ Phi : num -2.362 2.94 1.307 -0.381 -2.903 ...
> $ Theta : num NA 5.3 -1.63 -1.69 -2.52 ...
> $ Step.ID : num 1 2 3 4 5 6 7 8 9 10 ...
> $ Season : chr "winter" "winter" "winter" "winter" ...
> $ In.buffer : logi TRUE TRUE TRUE TRUE TRUE TRUE ...
> $ Dist_toBarrier: num 1702 4966 4928 2085 2095 ...
> $ Crossed : logi NA FALSE FALSE FALSE FALSE FALSE ...barrier contains the barrier data, with rows for every
segment and including columns for:
Complex locations (see as.complex()) for the start
and end of each barrier segment (Z1,
Z2)
Identifier for the each barrier segment
(barrier.id)
Identifier for each line (in sf terms, LINESTRING)
in the barrier (line_id) - will be set to 1 if there is
only 1 unique line
Any barrier-specific covariate values (km)
str(bor_format$barrier)
> 'data.frame': 4784 obs. of 5 variables:
> $ Z1 : cplx 472462+6771187i 472177+6771241i 472153+6771245i ...
> $ Z2 : cplx 472177+6771241i 472153+6771245i 472135+6771248i ...
> $ line_id : num 1 1 1 1 1 1 1 1 1 1 ...
> $ km : num 0.29 0.314 0.333 0.361 0.381 ...
> $ barrier.id: int 1 2 3 4 5 6 7 8 9 10 ...The print or summary functions can also be
used to get a quick summary of the permdata object,
including the final buffer size (Buffer.dist.m), number of
ID’s with locations in the buffer (N.IDs.in.Buffer), sample
rate and date range of the data (Sample.rate.hrs,
Date.range), number of total steps in the buffer
(N.steps.in.buffer), the number of steps that crossed the
barrier (N.steps.to.Cross), and the mean steplength for
steps that crossed the barrier
(Mean.crossing.steplength.m).
Null steps for each actual, observed movement step are simulated by sampling observed step lengths and angles from the data.
The buildCrossingTable function will default to creating
null steps by:
sampling step lengths for all steps that cross the barrier
sampling all turning angles for steps within the defined buffer
distance of the barrier, using n.null as the sample
size
Notably, if there are < 30 steps that cross the barrier, a Weibull
distribution will be fit to the steps that cross, using the mean
crossing steplength as a starting scale value (with the
fitdistr function from the MASS package); the
rweibull function will then be used to sample
n.null number of steplengths from this distribution.
Users can also optionally provide their own vectors of step lengths
and angles (with the arguments: step.lengths,
turning.angles, absolute.angles) to generate
null steps. Users should ensure that these accurately represent steps
with the “potential to cross” a given barrier, however, as
fitPermeability is sensitive to the null distribution.
The crossing table can be created for all individuals at once using
the buildCrossingTable function.
Main input arguments are the:
permdata list object with track and
barrier dataframes (output from the
prepPermeability function)
Number of null steps to be simulated for each actual step
(n.null, default 60)
Note that this function works by determining all the stepwise crossings for each animal movement step and barrier segment within the buffer and at the desired sample rate, so for large datasets, this function can take 10+ minutes to run.
For very large datasets, running in parallel is recommended to
significantly improve processing time. This option uses the
future.apply package.
The resulting crossing table, similar to the track
object in permdata as output from
prepPermeability, has rows for every movement step with
null crossings (e.g., includes every step that had the “potential to
cross”) and additional columns for:
barrier.id (barrier segment where the majority of
null steps crossed)
null.crossed and null.stayed (number of
null steps to cross versus not cross the barrier)
crossed (1 if the actual step crossed, 0 if
not)
barrier.crossed (barrier segment where the actual
step crossed - NA if it did not)
str(bor_cT)
> 'data.frame': 2473 obs. of 10 variables:
> $ ID : Factor w/ 10 levels "ID_1","ID_11",..: 1 1 1 1 1 1 1 1 1 1 ...
> $ Step.ID : num 2 3 4 5 6 7 8 9 10 11 ...
> $ barrier.id : num 984 982 993 989 989 989 989 993 989 989 ...
> $ Z.start : cplx 327542+6784926i 327254+6784985i 328011+6787781i ...
> $ Z.end : cplx 327254+6784985i 328011+6787781i 328043+6787769i ...
> $ Time : POSIXct, format: "2015-02-21" "2015-02-21" ...
> $ null.crossed : num 7 10 15 25 16 23 22 22 23 25 ...
> $ null.stayed : num 53 50 45 35 44 37 38 38 37 35 ...
> $ crossed : num 0 0 0 0 0 0 0 0 0 0 ...
> $ barrier.crossed: num NA NA NA NA NA NA NA NA NA NA ...In order to fit permeability models with covariate(s),
the crossing table needs to be annotated with any of these
variables.
The crossing table has rows corresponding to the interaction of each
movement step (identified with Step.ID) and each barrier
segment (identified with barrier.id). It additionally
contains a Time column (in POSIXct format) that contains
the timestamp for the first location in the movement step, which can be
used to match with time-varying covariates.
For the caribou data, covariates of interest include season (categorical), traffic volume (numeric), day of year (numeric), and kilometers along the highway (numeric).
Season can be added to the crossing table by simply merging the
crossing table with its permdata object, using the matching
ID and Step.ID columns.
track_covar <- bor_format$track[,c("ID", "Step.ID","Season")]
bor_cT <- merge(bor_cT, track_covar, by=c("ID", "Step.ID"),
all.x=TRUE)Kilometer markers along the highway is also available in the
permdata object, specifically in the barrier
dataframe (km). It can be added to the crossing table using
the matching barrier.id values.
barrier_covar <- bor_format$barrier[,c("barrier.id", "km")]
bor_cT <- merge(bor_cT, barrier_covar, by=c("barrier.id"),
all.x=TRUE)Traffic volume is contained in a separate dataframe, with daily
traffic volumes recorded for each unique date. Date can be
extracted the from the Time column in the crossing table
and used to merge with the traffic data matching Date
column.
bor_cT$Date <- as.Date(bor_cT$Time)
bor_cT <- merge(bor_cT, traffic_data, by = "Date", all.x = TRUE) |>
unique()Finally, day of year can be extracted using the
lubridate package and its yday function on the
Time column in the crossing table.
str(bor_cT)
> 'data.frame': 2473 obs. of 15 variables:
> $ Date : Date, format: "2015-02-21" "2015-02-21" ...
> $ barrier.id : num 984 982 989 993 989 989 989 993 984 989 ...
> $ ID : Factor w/ 10 levels "ID_1","ID_11",..: 1 1 1 1 1 1 1 1 1 1 ...
> $ Step.ID : num 2 3 6 4 5 7 8 9 12 10 ...
> $ Z.start : cplx 327542+6784926i 327254+6784985i 328004+6787759i ...
> $ Z.end : cplx 327254+6784985i 328011+6787781i 328463+6788763i ...
> $ Time : POSIXct, format: "2015-02-21" "2015-02-21" ...
> $ null.crossed : num 7 10 16 15 25 23 22 22 23 23 ...
> $ null.stayed : num 53 50 44 45 35 37 38 38 37 37 ...
> $ crossed : num 0 0 0 0 0 0 0 0 0 0 ...
> $ barrier.crossed: num NA NA NA NA NA NA NA NA NA NA ...
> $ Season : chr "winter" "winter" "winter" "winter" ...
> $ km : num 157 157 159 159 159 ...
> $ traffic_volume : num 53 53 55 55 55 59 59 59 55 55 ...
> $ doy : num 52 52 53 53 53 54 54 54 55 55 ...The plotCrossings function can be used with the
permdata object to visualize crossings, with the default
being a base R plot for all individuals. Setting
individual.plots=TRUE will output plots for tracks
individually.
More advanced options for plotCrossings include an
interactive map (interactive.map=TRUE, using the
mapview package) or a static map
(static.map=TRUE, ggplot2 package, with the
option to include a basemap with add.basemap=TRUE using the
ggspatial package). Static maps can be easily exported as
images.
Importantly, to use the interactive and static options for
plotCrossings, the user also needs to provide the original
data and the sf version of the barrier, for spatial
mapping.
plotCrossings(permdata = bor_format,
original.data = boreal_caribou,
barrier.sf = hwy1_west,
interactive.map = TRUE)crossing_map <- plotCrossings(permdata = bor_format,
original_data = boreal_caribou,
barrier.sf = hwy1_west,
static.map = TRUE,
add.basemap = TRUE)
crossing_mapThe null model can be fit to the crossing table data to estimate \(\kappa\) for the barrier of interest, similar to the simulation example.
>
> Estimate of kappa (with 95% confidence intervals):
> kappa.hat ci.low ci.high
> 1 0.04057158 0.02762201 0.05959209
Covariate models can include additive and/or interactive fixed effects and/or optionally, a splined covariate.
First, a model with traffic volume:
> Permeability fit call: ~traffic_volume
>
> estimate se z_score p_value
> (Intercept) -3.2341552 0.2023819 -15.980457 0.0000000
> traffic_volume -0.3093507 0.2012120 -1.537437 0.1241864
>
> logLik: -131.8029
> AIC: 267.6058
Second, a model with season as a categorical variable:
> Permeability fit call: ~Season
>
> estimate se z_score p_value
> (Intercept) -3.27469848 0.2094952 -15.63137717 0.0000000
> Seasonfall -0.01309219 0.2570990 -0.05092276 0.9593871
> Seasonpre-calving 0.19901172 0.1688073 1.17892832 0.2384267
> Seasonsummer -0.33995220 0.3118657 -1.09005957 0.2756869
> Seasonwinter -0.07383229 0.2810450 -0.26270630 0.7927770
>
> logLik: -130.4724
> AIC: 270.9447
Models with interactive and additive effects of season can then be specified:
> Permeability fit call: ~Season + traffic_volume
>
> estimate se z_score p_value
> (Intercept) -3.24351343 0.2079164 -15.60008260 0.0000000
> Seasonfall 0.04886779 0.2920432 0.16733070 0.8671098
> Seasonpre-calving 0.25741757 0.2004034 1.28449715 0.1989681
> Seasonsummer -0.28763592 0.3272580 -0.87892713 0.3794408
> Seasonwinter 0.03767124 0.3609408 0.10436958 0.9168761
> traffic_volume -0.00920928 0.3202563 -0.02875597 0.9770592
>
> logLik: -130.5901
> AIC: 273.1802
> Permeability fit call: ~Season * traffic_volume
>
> estimate se z_score p_value
> (Intercept) -3.3738979 0.237806 -14.1876067 0.0000000
> Seasonfall 2.8282825 2.475522 1.1424994 0.2532465
> Seasonpre-calving 2.4536024 1.751745 1.4006615 0.1613153
> Seasonsummer 1.1038139 3.059787 0.3607487 0.7182874
> Seasonwinter 0.4360159 2.996915 0.1454882 0.8843254
> traffic_volume 1.1409601 1.165786 0.9787049 0.3277258
> Seasonfall:traffic_volume -2.2562697 1.986646 -1.1357182 0.2560745
> Seasonpre-calving:traffic_volume -1.7571677 1.428196 -1.2303403 0.2185697
> Seasonsummer:traffic_volume -1.3729088 3.082937 -0.4453250 0.6560849
> Seasonwinter:traffic_volume 0.2784876 2.130988 0.1306847 0.8960247
>
> logLik: -130.6696
> AIC: 281.3392
Models with a spline of a continuous numeric covariate can also be
specified, e.g., day of year or kilometers along the highway. The
fitPermeability function uses a wrapper function with
mgcv functions for fitting smoothers and has support for
mgcv arguments (e.g., number of knots, smoother type, etc),
using the s function.
Here, we specify number of knots (k) of 4 and cyclic
cubic regression spline type (bs = "cc") for day of year as
a covariate.
> Permeability fit call: ~s(doy, k = 4, bs = "cc")
>
> Linear coefficients:
>
> estimate se z_score p_value
> Intercept -3.237283 0.2026308 -15.97627 0
>
> Approximate significance of smooth terms:
>
> k LRT p_value
> s(doy) 3 2.363091 0.5005422
>
> logLik: -131.839
> AIC: 271.678
Results include the number of spline parameters (k, also
the number of basis functions) and results of a Likelihood Ratio Test
(including test statistic and p-value) comparing the model with the
splined covariate to a model without it (in this case, just a null model
of permeability). A p-value < 0.05 indicates that the model with the
splined covariate is significantly improved over one without it.
Spatial variability in permeability along a highway can be determined using a splined covariate of kilometers along the highway:
> Permeability fit call: ~s(km, k = 6, bs = "cs")
>
> Linear coefficients:
>
> estimate se z_score p_value
> Intercept -3.515376 0.5342036 -6.580593 4.685741e-11
>
> Approximate significance of smooth terms:
>
> k LRT p_value
> s(km) 6 12.73855 0.04738078
>
> logLik: -126.6513
> AIC: 267.3025
Additive effects can also be used with a splined covariate:
doy_traffic_model <- fitPermeability(~ traffic_volume + s(doy, k = 4, bs = "cc"),
data = bor_cT)
doy_traffic_model> Permeability fit call: ~traffic_volume + s(doy, k = 4, bs = "cc")
>
> Linear coefficients:
>
> estimate se z_score p_value
> (Intercept) -3.2416058 0.2038093 -15.9050925 0.0000000
> traffic_volume -0.2778513 0.3619453 -0.7676609 0.4426887
>
> Approximate significance of smooth terms:
>
> k LRT p_value
> s(doy) 3 0.5589067 0.9057718
>
> logLik: -131.5234
> AIC: 273.0469
Models can be easily compared using the comparePermFits
function from permeability, which easily extracts and sorts
AIC values for a set of candidate models.
Users can either provide a list of formulas or model objects (as
output from fitPermeability) to the function.
bor_models <- list(Null = null_model,
Traffic = traffic_model,
Season = season_model,
SeasonandTraffic = additive_model,
SeasonxTraffic = interactive_model,
SplineDoy = spline_doy_model,
SplineKm = spline_km_model,
SplineDoyandTraffic = doy_traffic_model)
compare_models <- comparePermFits(models = bor_models)Or:
bor_formulas <- list(Null = ~1,
Traffic = ~traffic_volume,
Season = ~Season,
SeasonandTraffic = ~Season + traffic_volume,
SeasonxTraffic = ~Season*traffic_volume,
SplineDoy = ~ s(doy, k = 4, bs = "cc"),
SplineKm = ~ s(km, k = 6, bs = "cs"),
SplineDoyandTraffic = ~ traffic_volume + s(doy, k = 4, bs = "cc"))
compare_models <- comparePermFits(bor_formulas, data = bor_cT)
compare_models> Model k logLik AIC dAIC
> 1 SplineKm 7 -126.6513 267.3025 0.0000000
> 2 Traffic 2 -131.8029 267.6058 0.3032416
> 3 Null 1 -133.0205 268.0411 0.7385545
> 4 Season 5 -130.4724 270.9447 3.6422041
> 5 SplineDoy 4 -131.8390 271.6780 4.3754633
> 6 SplineDoyandTraffic 5 -131.5234 273.0469 5.7443348
> 7 SeasonandTraffic 6 -130.5901 273.1802 5.8776728
> 8 SeasonxTraffic 10 -130.6696 281.3392 14.0366378
The output includes the number of parameters (k), log
Likelihood, AIC, and delta AIC (dAIC) values for each
model, ordered by lowest AIC to highest.
Here, there is less than 1 AIC unit difference between a null model of permeability, a model with a spline of km along the highway, and a model with traffic volume as a covariate, suggesting that the spline of km model, despite having the lowest AIC, is not significantly improved over the other two possible models.
The predict function (with se.fit = TRUE
for 95% confidence intervals) can be used with either new or existing
data to get predicted \(\kappa\) values
as a function of covariates.
We will make predictions for the top-ranking covariate model, with
traffic_volume.
predict_kappa <- data.frame(traffic_volume = traffic_model$data$traffic_volume)
prediction <- predict(traffic_model, se.fit = TRUE)
predict_kappa <- cbind(predict_kappa, prediction)The user can also provide new data with a range of possible covariate values to predict on:
predict_kappa <- data.frame(traffic_volume =
seq(0, 100, by = 1))
prediction <- predict(traffic_model, new.data = predict_kappa, se.fit = TRUE)
predict_kappa <- cbind(predict_kappa, prediction)Predictions (with 95% confidence intervals) can be plotted to visualize the relationship between \(\kappa\) and the covariate (daily traffic volume):
plot(predict_kappa$traffic_volume, predict_kappa$kappa.hat, type = "l",
xlab = "Daily Traffic Volume",
ylab = "Estimated Kappa Values",
col ="orange", lwd = 2, ylim = c(0, 0.15))
lines(predict_kappa$traffic_volume, predict_kappa$ci.high, col = "grey")
lines(predict_kappa$traffic_volume, predict_kappa$ci.low, col = "grey")The permeability package also supports boostrapping to
get standard errors on estimates, in the case of the Hessian being
unable to get standard errors. The fitPermeability and
predict functions will give warnings if standard errors are
unable to be estimated with the Hessian and will suggest that users use
bootstrap = TRUE with the predict function to
get standard errors and confidence intervals on estimates.
Bootstrapping can take awhile, so it’s also suggested that users run
the predict function with parallel = TRUE
(using the future.apply package) to speed up run time.
Predictions from fitted models can also be mapped to sf
barrier segments, using the mapPredictions function.
This function takes a dataframe of predictions, including a
barrier.id value for each value so that it can be mapped to
barrier segments, and the original sf barrier object used
as input to the prepPermeability function.
The function can also handle multiple linear features; for example,
if you had three highways of interest occurring in the same spatial
area, you could bind them together into one large sf
object, including a column Feature.ID to distinguish them.
The predictions dataframe should also have a column called
Feature.ID, to map individual linear feature predictions to
their own geometry.
Finally, this function has options for creating a static map without
a basemap (default), with a basemap (using ggspatial if
maptype is specified or basemaps with an ESRI
landscape map, by default), or to create an interactive map (with
mapview and interactive = TRUE). Users can
optionally add an sf POINT object of their movement data to
map movement tracks with the barrier(s) (using the
tracks.sf argument) or even add their own basemap (provided
to the basemap argument).
Here, we show an example predicting \(\kappa\) values as a function of
km along the highway, using the fitted
spline_km_model from above. First, a
predictions dataframe is created, storing km
and corresponding barrier.id values for each barrier
segment.
hwy1_west_z <- hwy1_west_xy[,1] + 1i*hwy1_west_xy[,2]
predictions <- cbind(hwy1_west_xy, hwy1_west_z) |>
data.frame() |>
group_by(line_id) |>
reframe(km = Mod(diff(hwy1_west_z))) |>
select(km) |>
mutate(km = cumsum(km)/1000) |>
mutate(barrier.id = 1:length(km)) |>
data.frame() |>
subset(barrier.id %in% # make predictions for barrier segments in crossing table (crossings possible)
seq(min(bor_cT$barrier.id), max(bor_cT$barrier.id, by = 1)))Predicted \(\kappa\) values are
added to the dataframe, using the predictions dataframe as
a new data input.
predictions$kappa.hat <- predict(spline_km_model, se.fit = FALSE,
new.data = predictions)[,1]
str(predictions)> 'data.frame': 3271 obs. of 3 variables:
> $ km : num 0.29 0.314 0.333 0.361 0.381 ...
> $ barrier.id: int 1 2 3 4 5 6 7 8 9 10 ...
> $ kappa.hat : num 0.142 0.142 0.142 0.142 0.142 ...
A basemap can be added with add.basemap = TRUE.
boreal_caribou_sf <- boreal_caribou |>
st_as_sf(coords = c("X","Y"), crs = st_crs(hwy1_west))
map <- mapPredictions(predictions, barrier.sf = hwy1_west,
add.basemap = TRUE,
tracks.sf = boreal_caribou_sf)
mapWe present a simple and fast approach for estimating permeability of linear barriers for large biotelemetry tracking datasets. Our approach provides similar flexibility to other statistical model packages and functions. Users can fit either univariate or multivariate models and can perform model selection to estimate barrier permeability, with the option to explore covariates specific to either the movement steps of the animal or the barrier.
This tool especially serves as a useful method for quantifying permeability with respect to a barrier itself, with broad applications for management of taxa interacting with anthropogenic infrastructure (or even natural features, such as rivers).
If you use our package in publications, please be sure to cite it as:
Note: will be available upon publication on CRAN
Report any issues or bugs to the package developers:
(Anonymous)
Boreal caribou tracking data used in this vignette was provided by the GNWT and made possible by the support of Indigenous Government Organizations and management authorities for boreal caribou, including the Sambaa K’e Dene Band, Fort Simpson Métis Local, Łı́ı́dlıı Kų́ę́ Kue First Nation, Tthets’ék’ehdélı First Nation, Pehdzeh Ki First Nation, Nahɂą Dehé Dene Band, Acho Dene Koe Band, Ka’a’gee Tu First Nation, Deninu K’ue First Nation, NWT Metis Nation, Hay River Métis Council, Fort Resolution Métis Council, K’atlodeeche First Nation, West Point First Nation, Deh Gáh Got’ıę First Nation, and Fort Providence Métis Council.